/*
 * File:   FEMCalcBase.h
 *
 *      Author: Jerome RD Soine and Christoph A Brand
 *      Institution: Schwarz goup, Institute for Theoretical Physics, Heidelberg University, Germany
 */

#ifndef FEMCALCBASE_H
#define	FEMCALCBASE_H

#include <lac/vector.h>
#include <lac/full_matrix.h>
#include <lac/sparse_matrix.h>
#include <lac/solver_cg.h>
#include <lac/solver_richardson.h>
#include <lac/precondition.h>
#include <lac/constraint_matrix.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/grid_refinement.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/tria_boundary_lib.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <fe/fe_values.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
#include <numerics/error_estimator.h>
#include <grid/grid_in.h>
#include <grid/grid_out.h>
#include <grid/grid_tools.h>
#include <base/timer.h>
#include <numerics/fe_field_function.h>
#include <base/function.h>
#include <base/geometry_info.h>
				
#include <fe/fe_system.h>
#include <fe/fe_q.h>

				 
#include <fstream>
#include <iostream>

#include "../../../core/ForceTransmission.h"

#include "../../../core/SubstrateModel.h"

using namespace dealii;

/*!
 *  Super class for FEM calculations.
 *
 *  This class provides basic variables and functions for finite element calculations. It uses the library deal.II.
 *  Own continuum model implementations shall be derived from this class.
 */

class FEMCalcBase : public SubstrateModel
{
public:
	/*!
	 * Class constructor with empty arguments
	 */
	FEMCalcBase();
	/*!
	 * Copy constructor.
	 */
	FEMCalcBase(FEMCalcBase const &);
	/*!
	 * Class destructor
	 */
	virtual ~FEMCalcBase();
	/*!
	 * Generates a rectangular grid based on two corner points p1. Steps of local adaptive mesh refinements
	 * are based on a given ForceTransmission object in order to achieve high local refinement where force application
	 * occurs.
	 */
	void generate_rectangular_mesh(Point<3>,Point<3>,vector<ForceTransmission*> force_transmissions,unsigned int);
	/*!
	* Import mesh in GMSH format.
	*/
	void mesh_import(std::string);
	/*!
	* Import mesh in GMSH format and apply refinement history
	*/
	void mesh_import_with_history(std::string,std::string);
	/*!
	* Export current mesh as GMSH file.
	*/
	void mesh_export(std::string);
	/*!
	 * Writes displacement solution as VTK formated file.
	 */
	void output(std::string);
	/*!
	 * Writes force distribution as VTK formated file.
	 */
	void output_force(std::string);
	/*!
	 * Write both displacement and force into a given folder.
	 */
	void write_output (string filename) {output(filename+"_disp.vtk"); output_force(filename+"_force.vtk");};
	/*!
	 * Returns the number of DoFs associated with the displacement field.
	 */
	unsigned int n_dofs();
	/*!
	 * Returns the number of DoFs associated with the applied force field.
	 */
	unsigned int n_force_dofs();
	/*!
	* Returns the number of vertices.
	*/
	unsigned int n_vertices();
	/*!
	 * Apply global refinement
	 */

	void refine_globaly();
  void refine_surface();
  void refine_adaptively();
  void set_refine_parameter(double);
  
  void reinitialize_values(Vector<double>);
  
  void point_evaluation(std::vector<Point<3> > &points, std::vector<Vector<double> > &value_list);
  void point_evaluation(std::vector<Point<3> > &points, std::vector<Vector<double> > &value_list, std::vector<DoFHandler<3>::active_cell_iterator>& cell_iterators);
  void point_evaluation(Point<3>& p, Vector<double>& v);
  
  void save_solution_vector(std::string);
  void save_force_vector(std::string);
  
  DoFHandler<3>::active_cell_iterator find_cell_around_point(Point<3> p);
  
  Vector<double> get_solution();
  void set_solution(Vector<double>);
  
  void set_force_solution(Vector<double>);

  void write_force_field(std::string);
  
  Vector<double> get_force_vector();
  
  Vector<double> get_force_solution_vector();

  double force_sum();
  
 // virtual void run () = 0;
  virtual void run (vector<ForceTransmission*> force_trans) = 0;
  
  vector<double> get_displacement_at (vector<double> point);

  vector<vector<double> > get_displacement_list_at (vector<vector<double> > points);

  Triangulation<3> get_triangulation();

  double get_total_traction();

  void project_traction(vector<ForceTransmission*> force_trans);

protected:
  void reinit_force();
  virtual void setup_system(){};
  virtual void assemble_system() {};
  virtual void assemble_rhs() {};
  virtual void apply_boundary_conditions() {};
  virtual void solve() {};
  
  void refine_based_on_traction_shapes(vector<ForceTransmission*> force_transmission);

  vector<double> point_to_std_vector(dealii::Point<3> &point);



  virtual void move_mesh();

  Triangulation<3>* triangulation;
  
  DoFHandler<3>* dof_handler;
  DoFHandler<3>* dof_handler_force;

  FESystem<3>* fe;
  FESystem<3>* fe_force;
  
  ConstraintMatrix     hanging_node_constraints;
  SparsityPattern      sparsity_pattern;
  SparseMatrix<double> system_matrix;
  SparseMatrix<double> recover_matrix;

  Vector<double>       solution;
  Vector<double>       system_rhs;

  std::vector<DoFHandler<3>::active_cell_iterator> cell_iterator_list;
  bool first_calculation;

  bool first_force_projection;
  vector<DoFHandler<3>::active_face_iterator> traction_faces; // Save face iterator associated with forces
  vector<vector<unsigned int> > face_to_force_indices; // Save forcetrans indices associated with faces above

  Vector<double>       force;
  std::vector<unsigned int> force_indices;
  Vector<double>       force_solution;
  
  double total_area;
  double cell_area;


  double ref_parameter;
  
  static const unsigned int boundary_index = 5;
  
  unsigned int n_forces;
  
};
  


#endif
